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Abstract. In two-dimensional reaction-diffusion systems, local curvature 
perturbations in the shape of traveling waves are typically damped out and 
disappear in the course of time. If, however, the inhibitor diffuses much faster 
than the activator, transversal instabilities can arise, leading from flat to folded, 
spatio-temporally modulated wave shapes and to spreading spiral turbulence. 
For experimentally relevant parameter values, the photosensitive Belousov- 
Zhabotinsky reaction (PBZR) does not exhibit transversal wave instabilities. 
Here, we propose a mechanism to artificially induce these instabilities via a 
wave shape dependent spatio-temporal feedback loop, and study the emerging 
wave patterns. In numerical simulations with the modified Oregonator model 
for the PBZR using experimentally realistic parameter values we demonstrate 
the feasibility of this control scheme. Conversely, in a piecewise-linear version 
of the FitzHugh-Nagumo model transversal instabilities and spiral turbulence in 
the uncontrolled system are shown to be suppressed in the presence of control, 
thereby stabilising flat wave propagation. 
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1. Introduction 

A large variety of pattern forming processes can be understood in terms of the 
advancement of an interface between two or more spatial domains. An interface that 
becomes unstable to diffusion possibly causes intricate spatio-temporal dynamics. Well 
known examples include the Mullins-Sekerka instability during crystal growth and 
formation of snow flakes [1, 2], and the Saffman-Taylor instability leading to viscous 
fingering in multiphase flow and porous media [3, 4, 5]. Other phenomena affected by 
interfacial instabilities are flame fronts [6, 7], Marangoni convection [8], and growing 
cell monolayers [9]. 

Traveling plane waves in excitable media exhibit interfacial instabilities as well. Here, 
an effective interface separates the excited state from the excitable rest state. A 
straight isoconcentration line of a two-dimensional flat wave can suffer an instability 
leading to stationary or time dependent modulations orthogonal to the propagation 
direction. Further away from the instability threshold, rotating wave segments and 
spreading spiral turbulence are observed [10, 11]. For standard activator-inhibitor 
kinetics, these so-called transversal or lateral wave instabilities typically occur if the 
inhibitor diffuses much faster than the activator. This result was analytically predicted 
first by Kuramoto for piecewise-linear reaction kinetics [12, 13]. Later, it was confirmed 
numerically by Horvath et al. for autocatalytic reaction-diffusion fronts with cubic 
reaction kinetics [14] as well as in experiments with the iodate-arsenous acid reaction 
[15] and the acid-catalyzed chlorite-tetrathionate reaction [16]. 

The experimental workhorse of chemical pattern formation, the Belousov-Zhabotinsky 
(BZ) reaction, does typically not display transversal wave instabilities. Dispersing 
the reagents of the BZ reaction in nanodroplets of a water-in-oil microemulsion 
allows to increase the inhibitor diffusivity considerably [17] and leads, for example, to 
segmented spiral waves as reported by Vanag and Epstein [18]. Even in the presence 
of an electrical field applied to enhance transversal instabilities in cubic autocatalytic 
reaction-diffusion fronts, the inhibitor diffusion coefficient is always required to be 
sufficiently larger than that of the activator [19, 20, 21]. 

Because of the possibility to apply spatio-temporal external forcing or feedback- 
mediated control loops by exploiting the dependence of the local excitation threshold 
on the intensity of applied illumination, the photosensitive variant of the BZ reaction 
has been widely used as a paradigm of an experimentally well controllable RD 
system. So far, unstable wave propagation has been stabilised by global feedback 
[22]. Two feedback loops were used to stabilise unstable wave segments and to 
guide their propagation along pre-given trajectories [23]. Also, spiral wave drift in 
response to resonant external forcing and various feedback-mediated control loops 
has been extensively studied experimentally in PBZR systems, compare for example 
[24, 25, 26, 27, 28]. 

In this paper, we design a curvature-dependent spatio-temporal feedback loop in 
order to destabilise a stable propagating planar reaction-diffusion wave by inducing 
transversal instabilities. In numerical simulations with the modified Oregonator model 
for the PBZR, we study the wave patterns emerging beyond the instability threshold, 
and demonstrate the capability to actively select wave patterns by modifying 
feedback parameters accessible to the experimenter. Conversely, under conditions 
where planar wave propagation fails due to transversal instabilities, using the same 
feedback mechanism we suppress ongoing breakup and segmentation of waves, thereby 
stabilising unstable propagating planar waves. 
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2. Theory 
2.1. Models 


The first model we investigate is the three component modified Oregonator model [29] 

(1) 
( 2 ) 
( 3 ) 
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Here, the parameters e and e characterise the time scales for the dynamics of 
the activator u and inhibitor re, respectively, and the stoichiometric parameters 
q and / depend on the temperature and chemical composition. All parameter 
values used in numerical simulations are listed in table Al in Appendix A. The 
modified Oregonator model describes the light-sensitive Belousov-Zhabotinsky (BZ) 
reaction. In experiments, the catalyst v can be immobilised in a gel and therefore 
the corresponding diffusion coefficient is set to zero. The activator u and inhibitor 
w diffuse with diffusion coefficients Du and whose ratio for typical BZ recipes is 
approximately Dyj/Du ~ 1.2. This value is too low to support transversal instabilities 
such that plane waves are stable for typical BZ recipes. The parameter ^ in equation 
(3) is proportional to the applied light intensity and measures the local excitation 
threshold. In experiments, spatio-temporal modulations of ^ can be applied to control 
wave propagation in the BZ reaction [23, 27, 28, 22]. 

Because the modified Oregonator model does not exhibit transversal instabilities in 
the parameter regime relevant for experiments, we investigate a second model. The 
piecewise-linear caricature of the FitzHugh-Nagumo (FHN) model [30, 11] received 
some attention in the context of transversal instabilities [11]. It is a two component 
model of standard activator-inhibitor type. 


du ^ ^, 

— = Au + F{u,v), 
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with u being the activator and v the inhibitor. The reaction kinetics are a piecewise- 
linear caricature of the FHN model 


where 


F{u,v) = f{u) -V, 

G {u, v) = kgU — V, 

{ —kiu^ u < 5, 

kf {u — a), 6 < u < 1 — 6, 

k 2 {l — u), 1 — 6 < u. 


(6) 

( 7 ) 


(8) 


The parameters ki and k 2 are chosen such that / (u) is continuous ai u = 6 and 
u = 1 — S, which leads to 


ki = ^(a-S), 


^2 = h (1 - (5 - a). 


( 9 ) 
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Figure 1. Slightly beyond the onset of transversal instabilities, an initially plane 
wave in the piecewise linear FitzHugh-Nagumo model develops a stationary fold. 
From left to right, top to bottom: time sequence of snapshots of a video which 
can be found in the Supplemental Material [31] for (a) t = 100, (b) t = 320, 
(c) t = 370, (d) t = 420, (e) t = 470, and (f) t = 570 for the parameters 
a = 2.1, e = 0.1425, a = 0.1. Computations were carried out with a constant 
time step At = 0.00001 and space step Ax = 0.15. The domain size is 160 x 90. 


The remaining parameters for the function / are chosen in such a way that / resembles 
the cubic shape of the FHN activator nullcline. All parameter values used in numerical 
simulations are listed in table A2 in Appendix A. The parameter a is a measure for 
the excitation threshold and used as the feedback parameter. 

For numerical simulations, we assume an elongated two-dimensional channel of width 
L in the ^-direction with waves propagating in the x-direction. The boundary 
conditions in the x-direction are periodic while we assume periodic or Neumann 
boundary conditions in the ^-direction. For both models, we use a box-like initial 
condition of width b for the vector of components u. 


u {x, y, to) = ©Box {{x - (p (y, to)) /b) (u, 


uo) Uo, 


( 10 ) 


Lmax 


where Umax is the initial height of the pulse and Uq is the stationary point of the 
reaction kinetics. The box function is defined as 



( 11 ) 


The initial shape of the wave is given by 



( 12 ) 


where A denotes the amplitude of deviation of the shape from a plane wave and d is 
an offset. For numerical simulations in two spatial dimensions, we use Euler forward 
for the time evolution and a five point stencil for the Laplacian. 

A phase diagram for the occurrence of transversal instabilities in the e-cr-parameter 
plane of the piecewise-linear FHN model was presented by Zykov et al. in [11]. In¬ 
creasing the inhibitor diffusion coefficient a crosses the threshold for transversal in¬ 
stabilities. Shortly beyond the onset of transversal instabilities, a plane wave develops 
a fold which is stationary in a comoving frame of reference, see figure 1 for a time 
sequence of snapshots and the supplemental material [31] for a movie. Further away 
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from the instability threshold, a plane wave breaks into segments which undergo self- 
sustained rotatory motion accompanied by permanent merging and annihilation of 
segments. This regime is also known as spreading spiral turbulence [11], see figure 2 
for a time sequence of snapshots. 


2.2. Evolution equations for isoeoncentration lines 

Theoretically, the onset of transversal instabilities can be understood with the help of 
the linear eikonal equation 

Cn {S,t) = C - UK{s,t) , (13) 

an evolution equation for a two-dimensional curve 7 (s,t) = (s, t), (s, t))^ 

The 


representing an isoconcentration line parametrised by the curve parameter s. 
linear eikonal equation relates the normal velocity (n is the normal vector of y) 

Cn {s,t) =n{s,t) ■ dtnf{s,t) 
along y linearly to its curvature, 

dsJx {s, t) d^'Yy (s, t) - ds'jy (s, t) 5^7^ (s, t) 


n{s^t) = 


{S,t)f + {ds'jy 


3/2 


(14) 


(15) 


The curvature is conventionally assumed to be positive for convex isoconcentration 
lines, i.e., an isoconcentration lines with a protrusion in the propagation direction. The 
constant c corresponds to the pulse velocity of a one-dimensional solitary wave and n 
is the curvature coefficient. A rigorous derivation of the eikonal equation (13) from 
the reaction-diffusion system identifies the constant n in terms of the one-dimensional 
pulse profile, its response function and the matrix of diffusion coefficients, see [32] for 
details. For a plane wave, any isoeoncentration level is a straight line and therefore 
its curvature vanishes, k (s, t) = 0, everywhere along y. The stability of a plane wave 
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Figure 2. Segmentation of initially plane waves and spreading spiral turbulence 
occurs deep in the regime of transversal instabilities in the piecewise linear 
FitzHugh-Nagumo model. Parameter values as in figure 1 except e = 0.1575. 
The segments undergo self-sustained rotatory motion and nucleate new waves. 
From left to right, top to bottom: snapshots of a video [31] for (a) t = 5, (b) 
t = 301, (c) t = 500, (d) t = 672, (e) t = 700, and (f) t = 891. The domain size 
is 160 X 90. 
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is determined by the sign of the curvature coefficient v. As long as z/ > 0 , any point 
along the isoconcentration line of a convex bulge moves slower than a plane wave, 
while points of a concave dent move faster than a plane wave, thereby smoothing out 
deviations from a plane wave. If < 0, a convex bulge moves faster than a plane wave, 
protruding the bulge even further and thereby leading to an ever increasing curvature: 
a transversal instability arises. Patterns arising for z^ < 0 cannot be described by the 
linear eikonal equation and terms depending nonlinearly on the curvature have to be 
taken into account which saturate the growth of an ever increasing curvature. At least 
two different nonlinear versions of equation (13) exist in the literature. Zykov et al. 
[33, 30, 34, 35] renormalised e and a in equation (5) to derive a renormalised one¬ 
dimensional velocity c depending on the curvature. Dierckx et al. [32] derive higher 
order nonlinear corrections in the curvature by a rigorous perturbation expansion with 
a small parameter proportional to the curvature, additionally generalising the eikonal 
equation to anisotropic media. 

Apart from nonlinear eikonal equations, which are difficult to solve numerically, 
patterns arising beyond the threshold of transversal instability can be described by 
the Kuramoto-Sivashinsky (KS) equation, 

dt4> iy,t) = c + I {dy(j) (y, t)f + vdl<j) (y, t) - Xd^cj) (y, t). (16) 

Equation (16) is an evolution equation for the x-component (j){y^t) of an 
isoconcentration line 7 parametrised in the form 7 (^, t) = (0 (^, t) ^y) . See [13] 
and [36] for a derivation of equation (16) from a general RD system. The case of 
Neumann boundary conditions in the ^-direction for the RD system imply that any 
iso concentration line of activator and inhibitor meets the domain boundary in a right 
angle. This corresponds to Neumann boundary conditions for 0, 

5^0 (0,t) =0, a^(/)(L,t) = 0. (17) 

Similarly, periodic boundary conditions in the ^-direction of the RD system carry 
over to periodic boundary conditions for (j). Equation (16) was originally proposed by 
Sivashinsky [ 6 ] in the study of turbulent flame propagation and adapted for reaction- 
diffusion systems by Kuramoto [37, 13]. The parameter A can be expressed in terms 
of a sum over all eigenfunctions of the linear stability operator arising through a 
linearisation of the one-dimensional RD system around the traveling wave solution 
[13]. To compute A, we use a method which avoids the virtually impossible numerical 
computation of all eigenfunctions, see [36] for details. The values of A and u for the 
modified Oregonator model with parameters as given in Appendix A are 

A = 0.68, z/ = 1.05. (18) 

The Kuramoto-Sivashinsky equation (16) allows a refined investigation of the onset of 
transversal instabilities. Eor a stability analysis of a plane wave in channel of width L 
with Neumann boundary conditions, we apply a perturbation expansion in 0 < e <C 1 
with an ansatz in form of a Eourier series, 

00 

4){y,t)=ct + e y] an exp {ujnt) cos (fff (19) 

n= —oo 

where the term ct corresponds to a plane wave solution to the RD system traveling in 
the x-direction. The dispersion relation follows as 
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Transversal instabilities occur only if uji > 0, i.e., v must be negative and the channel 
width must exceed 



Thus, in general, the onset of a transversal instability depends on the boundary 
conditions and can be suppressed in thin channels. It is a long-wavelength instability, 
i.e., the first mode which becomes unstable upon reaching the threshold is the mode 
with the longest possible wavelength. 

If z/ < 0, the fourth order term in the KS equation (16) counteracts the 
negative diffusion term and leads to a saturation of the growth of wavefront 
modulations.Starting at the threshold of instability, the solution to the KS 
equation(16) displays a fold with a minimum located at y = Lj^. Upon increasing 
L, this steady wave loses stability via a supercritical Hopf bifurcation [14] and the 
wave front starts to oscillate back and forth in a symmetrical fashion. Increasing L 
even further leads to a symmetry breaking bifurcation with asymmetrical oscillations 
followed by a period doubling cascade to fully developed spatio-temporal chaos. In 
this regime, the KS equation displays a strong dependence on the initial data, with 
small differences in the initial conditions leading to dramatically different future time 
evolution. This characteristic of the KS equation is also studied as an analogy for 
hydrodynamic turbulence [38]. 

As long as z/ > 0 , no instability can arise and the fourth order term can be safely 
neglected by setting A = 0. In this case, equation (16) simplifies to the nonlinear 
phase diffusion equation, which in turn can be transformed to the usual diffusion 
equation via the Cole-Hopf transform [ 12 ]. Therefore, equation (16) with A = 0 can 
be solved analytically for arbitrary initial and boundary conditions. 

To assess the accuracy of the Kuramoto-Sivashinsky equation (16) as an approximation 
for propagating reaction-diffusion waves, we compare the transition from an initially 
curved shape to a plane wave for z/ > 0 with numerical simulations of the underlying 
two-dimensional modified Oregonator model Eqs. (l)-(3). The isoconcentration line 
7 of the activator variable u is determined numerically as the set of points r = (x, y^ 
for which iz(r, t) = Uc = 0 . 2 . We compute the x-component of the wave’s centre of 
mass in a comoving frame as 

L 

(x) {t) = ^ J (j) {y, t)dy-(p (0, t). (22) 

0 

Figure 3 shows the time evolution of {x) (t) obtained from the Kuramoto-Sivashinsky 
equation (black dotted line) and nonlinear phase diffusion equation (blue solid line) 
and for the modified Oregonator model obtained by numerical simulations (red 
dashed line) for two different values of the amplitude A which characterises the 
initial deviation from a plane wave. As one would expect intuitively, the agreement 
between numerical simulations on the one hand and Kuramoto-Sivashinsky equation 
and nonlinear phase diffusion equation on the other hand becomes worse the larger is 
the initial amplitude A. For large times, i.e., when the curved isoconcentration line 
approaches a straight line, all results agree. The nonlinear phase diffusion equation 
and the Kuramoto-Sivashinsky equation practically yield the same result for all times, 
confirming the fact that the fourth order derivative in the Kuramoto-Sivashinsky 
equation can safely be neglected if the curvature coefficient is z/ > 0 . 
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Figure 3. Time evolution of the mean x-coordinate of an isoconcentration line 
in a comoving frame of reference. Comparison of Kuramoto-Sivashinsky equation 
(black dotted line), nonlinear phase-diffusion equation (blue line) and for the 
activator isoconcentration line with u (x, y, t) = 0.2 of the modified Oregonator 
obtained by numerical simulations (red dotted line). 



Figure 4. Velocity c of a one-dimensional solitary pulse over the parameter 4> 
proportional to applied light intensity for the modified Oregonator model. The 
result of numerical simulations (blue dots) can be well approximated by a linear 
least square fit (blue solid line). 


2.3. Curvature-dependent feedback control 

The feedback law proposed in this article requires that the velocity c of a one¬ 
dimensional wave depends sufficiently strongly on a parameter which can be controlled 
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Figure 5. Small modulations of the wave shape occur for weak feedback with 
an effective curvature coefficient of u = —0.75. Snapshots of a movie [31] for (a) 
t = 3.5, (b) t = 10.6, (c) t = 25.2, and (d) t = 69.3. Colours denote the value of the 
applied spatio-temporal illumination field ^ (a:, y, t) which is proportional to the 
curvature of the black isoconcentration line. Shown are clippings of size 21.5 X 30 
centred on the wave’s centre of mass while the computational domain is 110 X 30. 
This and the following computations of the modified Oregonator model were 
carried out with a time step size of At = 0.0001 and a space step size Ax = 0.05. 
Parameter values for the feedback law are = 0.018, ^max = 0.042. 

in experiments. For the modified Oregonator model, we use the parameter ^ 
proportional to the applied light intensity as the feedback parameter. A numerical 
computation of the dependence of c on ^ is shown in figure 4. With relatively good 
accuracy, the dependence can be assumed to be linear. 


c(4>) = Co + ci4>, (23) 

with parameters ci = — 90.191,co = 9.013 obtained from a least square fit. Solitary 
waves exist only in the excitable regime of 4> values indicated by the dashed lines in 
figure 4. For 4> < 0.045, the rest state is unstable and the medium becomes oscillatory. 
For 4> > 0.068, the solitary pulse profile becomes unstable and decays to the stable 
rest state. A successful feedback control is possible if <F is restricted to lie between 
these two values. 

We introduce a feedback law for 4> depending linearly on the curvature, 

4> (/^) = a fdn. (24) 

The parameters a and /3 are accessible to an experimenter. In general, these 
parameters can be adjusted with time to achieve a better performance of the control. 
Together with equation (23) and equation (24), the linear eikonal equation (13) 
becomes 

Cn = Cia + Co - nn, (25) 

with the effective curvature coefficient 

i> = 1 / — ci/3. 


( 26 ) 

















































Control of transversal instabilities in reaetion-diffusion systems 


10 


Depending on the sign of i>, the control will have very different effects. If a plane 
wave is stable with respect to transversal perturbations because z/ > 0 , we can excite 
transversal instabilities ifz> = z/ — ci/3 < 0 . Conversely, if z/ > 0 such that plane 
waves are unstable with respect to transversal modulations, patterns can be stabilised 
ifp = z/ — ci/3>0. An appropriate choice of the parameters a and /3 in the feedback 
law (24) allows to control transversal instabilities. 

To apply the feedback law (24) it is necessary to compute the curvature of a chosen 
iso concentration line of a chosen component with sufficient accuracy, which raises 
considerable difficulties. 


2.4- Computation of curvature by Level Set Methods 


The curvature k (s, t) of an isoconcentration line 7 ( 5 , t), equation (15), is proportional 
to the second derivative of the isoconcentration line with respect to the curve 
parameter s. Computations of isoconcentration lines from numerical simulations or 
experiments are affected by noise due to the discretised nature of the computed or 
measured concentration field iz, respectively. Numerical differentiation is an ill-posed 
mathematical operation and typically amplifies the noise. A variety of methods to 
compute the curvature n directly from a numerically determined isoconcentration line 
were tested and discarded due to insufficient performance. 

An indirect method which avoids the differentiation of an isoconcentration line is to 
compute the curvature field k as 

Viz (r, t) 


k (r, t) = V- 




(27) 


According to the formula of Bonnet [39], evaluating k at an isoconcentration line 
r = 7 (s, t) of u yields the curvature n of 7 , i.e.. 


K{s,t) (28) 

See Appendix B for a proof of Bonnet’s formula. Equation (27) involves the 
determination of the second derivative of u with respect to x and y. These expressions 
are readily available from the finite difference algorithm used to solve the RD system 
numerically. The problem is now that the concentration iz of a pulse solution typically 
varies very fast in a small spatial region while it is constant everywhere else, leading 
to an ill-defined denominator in equation (27). This difficulty can be addressed with 
the help of a level set method, which, however, is numerically quite expensive. 

Originally, level set methods were developed by Osher and Sethian to compute and 
track the motion of interfaces. These methods have since been successfully applied in 
such diverse areas of applications as computer graphics, medical image segmentation 
and crystal growth [40, 41, 42]. 

We introduce a second field variable y (r, r) which evolves in (virtual) time r according 
to the so-called reinitialisation equation [43, 42, 44] 


drX + sign (x°) (|Vx| - 1) = 0 

(29) 

( 1 if X > 0 

sign (x) = < 0 if X = 0 . 

1 —1 if X < 0 

(30) 


Equation (29) is solved with the initial condition 

X (r, 0) = x° (r) = u (r, t) - Uc, 


( 31 ) 
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Figure 6. Moderate feedback corresponding to an effective curvature coefficient 
of u = —1.059 leads to V-shaped patterns moving much faster than a plane wave. 
Snapshots of a movie [31] for (a) t = 2, (b) t = 7, (c) t = 62, and (d) t = 80. 
Shown are clippings of size 30 X 30 from a computational domain of size 110 x 30, 
and feedback parameter values are ^min = 0.018, ^max = 0.046. 

where Uc is the activator value along the isoconcentration line 7 for which we want 
to determine the curvature i.e., u (s^t) ,t) = Uc. Note that x( 7 (s,t),r) = 
(7 (s, t)) = 0 for all times r such that the position of the level set 7 is not changed 
by equation (29). However, equation (29) transforms the neighbourhood of x = 0 such 
that, after sufficiently many time steps r, 

lim |Vx(r,T)| = 1 . (32) 

r—)-oo 

The curvature of 7 , equation(15) can now readily be computed in terms of the 
Laplacian of x as 

K {s, t) = k{'y (s, t),t) = lim Ax (7 (s, t),T). (33) 

r—)-oo 

Numerically, the evolution of x ap to the final time r = 0.01 is sufficient to obtain a 
very accurate smooth result for the curvature of 7 . The reinitialisation equation (29) 
has to be solved at every real time step t. However, because the time evolution of the 
RD system is slow enough, we recompute the curvature k only every 200th time step 

t. 

3. Results 

3.1. Exeitation of transversal instabilities in the modified Oregonator model 

We study the possibility to excite transversal instabilities by curvature-dependent 
feedback in the modified Oregonator model. The feedback law (24) is realised via the 
parameter 4> proportional to the illumination in the BZ reaction. For the parameters 
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Figure 7. Segmentation of waves occurs under very strong feedback with an 
effective curvature coefficient of u = —2.337. Segments may break off, nucleate 
new waves and often start to rotate. Snapshots of a movie [31] for (a) t = 4, (b) 
t = 6, (c) t = 25, and (d) t = 50. The domain size is 110 X 30, and feedback 
parameters are 4>min = 0.006, 4>max = 0.05. 

of the feedback law we set a = ^max and P = — (^max — ^min) /^^norm such that the 
effective curvature coefficient is 

V — V C\P = n (^max ^min) • 

6 ^norm 

The values of ^h^ax and ^min can be chosen arbitrarily as long as ^min < ^max and 
both values lie in the regime of an excitable medium, see figure (4). The curvature k, is 
determined for the activator isoconcentration line 7 with (7 (s, t), t) = Uc = 0.2. An 
area of fixed size in front of and behind 7 is illuminated with the same value {n (s, t)), 
while within the remaining medium <l> attains its background value <l> = 4>o. Before 
the feedback is switched on at time ti = 0.4, the wave evolves uncontrolled. The value 
of /^norm = 1-2 is au estimate of the largest value which the curvature attains during 
the overall time evolution. For simplicity, we choose a constant value of /^norm, but in 
principle this value can be set to the maximum curvature every time the curvature is 
recomputed. 

Figure 5 shows wave patterns arising for weak feedback with an effective curvature 
coefficient i> = —0.75. The black solid lines denote the isoconcentration line 7 for 
the activator level Uc = 0.2. The rightmost line corresponds to the wave front while 
the trailing line corresponds to the wave back. The colours represent the value of 
the feedback parameter 4> and are proportional to the curvature of the wave front 
iso concentration line. An initially sinusoidal shape decays and a plane wave with 
transversal modulations of small wavelength develops. For the example presented 
here, the modulations are not stationary but travel along the isoconcentration line 
until they annihilate each other or at the Neumann boundaries. For even weaker 
feedback, the modulations do not travel such that the pattern is truly stationary in 
a comoving frame of reference. The overall velocity of the patterns is approximately 
the velocity c of the one-dimensional unperturbed traveling wave. Apart from the 
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^min 

Figure 8. Phase diagram for patterns in the modified Oregonator model under 
feedback control. Blue bullets • correspond to folded waves stationary in a 
comoving frame of reference, green triangles A stand for non-stationary wave 
front modulations, ochre squares ■ denotes segmentation of waves and spiral 
turbulence, a) refers to figure (7), b) to figure (6) and c) to figure (5). 


wave length of the modulations, this type of pattern appears similar to the patterns 
arising in the uncontrolled FHN model slightly beyond the threshold of transversal 
instabilities, see figure 1. 

Figure 6 displays the effects of moderate feedback with an effective curvature coefficient 
z> = —1.059. V-shaped patterns arise which travel much faster than a corresponding 
one-dimensional solitary pulse. In a frame of reference comoving with the centre 
of mass, the V-shaped patterns appear stationary apart from modulations traveling 
along the isoconcentration line. The V-patterns observed under feedback are long¬ 
time stable and do not decay or break up. A solitary V-pattern in an unbounded 
domain can be explained analytically as a solution to the linear and nonlinear eikonal 
equations [45, 46]. A V with opening angle a has a mean velocity V given by 


sin (a) ’ 

where c is the one-dimensional velocity. Because |sin((a)| < 1, all V-patterns are 
moving faster than a plane wave. Experimentally, these patterns were observed in 
homogeneous [47] and stratified [48] BZ media. 

Figure 7 shows the effect of strong feedback with an effective curvature coefficient 
L> = —2.337. Similar as for moderate feedback, V-shaped patterns appear. However, 
their shape is non-stationary but oscillating. The V-shape is segmented in an irregular 
and non-stationary way, with segments either merging again or breaking off and serving 
as the nucleation centre for new waves. These new waves propagate as segmented 
circles and occasionally start to rotate until they annihilate upon collision with other 
waves. Qualitatively, the segmentation and occurrence of rotating segments is similar 
to the spreading spiral turbulence observed for the uncontrolled FHN model deep in 
the regime of transversal instabilities, see figure 2. 

These results show that the proposed feedback law is not only able to excite transversal 
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Figure 9. Transversal instabilities are suppressed in thin channels. Shown is the 
stability boundary of a plane wave under feedback in a channel of width L with 
Neumann boundary conditions. Using feedback, the effective curvature coefficient 
u = u — ci^ is adjusted until a plane wave becomes unstable. Red line: theoretical 
prediction equation (21) obtained from the Kuramoto-Sivashinsky equation. Blue 
dots: result of numerical simulations of the modified Oregonator model. Blue line: 
least square fit to numerical results. 


instabilities but allows, to some extent, the selection of the patterns beyond the 
instability threshold by tuning the feedback parameters 4>niax and 4>niin accessible 
to an experimenter. We display a phase diagram with a classification of the observed 
patterns in the 4>niax — ^min plane in figure 8. Note that according to the KS equation 
(16), the observed patterns should only depend on the effective curvature coefficient 
given by equation (34). However, numerical simulations show that the type of 
patterns depends not only on the difference of 4>max and 4>min, but also display a slight 
dependence on their absolute values. This dependence is due to nonlinear corrections 
in the relation for the one-dimensional velocity c over 4> and higher order effects 
neglected by the KS equation (16). 

By adjusting the effective curvature coefficient z>, we are able to validate the predicted 
onset of transversal instabilities equation (21), 

L = (36) 

yJ-V 

and its dependence on the channel width L. We perform numerical simulations of 
the controlled Oregonator model in a channel with width L and Neumann boundary 
conditions in the ^-direction. Starting with a plane box-like initial condition equation 
(10) with noisy box width 6, we change the effective curvature coefficient z> until a plane 
wave becomes unstable, i.e., the curvature along the isoconcentration line is different 
from zero. Figure 9 shows that both numerical simulations and analytical prediction 
yield a linear relation between channel width L and \j \/—v over a large range of 
effective curvature coefficients z>. The slopes differ due to higher order corrections for 
the KS equation (16) and nonlinear corrections for the velocity c over 4>, equation 
(23), used for the feedback law. 

Beyond the onset of transversal instabilities, the emerging patterns can in principle 
be described by the KS equation (16). We compare the time evolution of the modified 
Oregonator model with the solution of the KS equation for an effective curvature 
coefficient of v = —0.02. Because the centre of mass velocity is incorrectly predicted 
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Figure 10. Time evolution of a wave’s isoconcentration line slightly beyond the 
threshold of transversal instability caused by feedback with an effective curvature 
coefficient 9 = —0.02. Compared are the results of two-dimensional numerical 
simulations of the modihed Oregonator model (dotted red line) with the solution 
of the Kuramoto-Sivashinsky equation (blue line). Due to the strong dependence 
on the initial conditions, the time evolutions are comparable only for a short time 
span, (a) corresponds to t = 0, (b) t = 6, (c) t = 12 and (d) t = 25. Comparison 
is made in the comoving frame of reference because the centre of mass velocities 
do not agree. 


by the KS equation, figure 10 shows a sequence of snapshots of isoconcentration lines in 
a frame of reference comoving with the centre of mass. Due to the strong dependence 
on the initial data, any initial agreement between the two curves is vanishing fast 
during the time evolution. 


3.2. Suppression of transversal instabilities 


The curvature-dependent feedback law (24) is used to suppress transversal instabilities 
occurring in the uncontrolled piecewise linear FHN model given by equations (4), (5). 
Here, the excitation threshold a is used as the feedback parameter. First, we linearly 
approximate the velocity - excitation threshold relation as 


c(a) = Co + cia. 


(37) 


with Co = 2.23 and ci = —8.75. Similar as shown in figure (4) for the modified 
Oregonator model, solitary pulses exist only for a certain range of a values. Second, 
the dependence of the excitation threshold a on the curvature k is chosen as 


rUniin T (d (t) /D, ^ 0, 

^Uniin <C 0. 


( 38 ) 
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The coefficient [3 (t) is adjusted in time such that the maximum value of a {k) along the 
iso concentration line does not exceed or undershoot the range of existence of solitary 
pulses. Every 100 time steps, we determine the maximum curvature /^max {t) along 
the isoconcentration line and set /^max to this value, 


/3(i) 


Ojmarx Ojx\ 


it) 


(39) 


The background value of a is set to uq = 0.1 everywhere before the feedback control 
is switched on at time t = ti = 215, and outside the region affected by the feedback 
control. Figure 11 displays the suppression of a transversal instability slightly beyond 
the threshold. For the same parameter values as in figure 1, the initially sinusoidally 
shaped wave relaxes back to a plane wave and no fold appears, see also the video in the 
supplemental material [31]. Patterns deep in the regime of transversal instabilities are 





u 


Figure 11. Curvature-dependent feedback control stabilises an unstable plane 
wave. For the same parameters slightly beyond the threshold of transversal 
instabilities as the corresponding uncontrolled time evolution in figure 1, the 
initially sinusoidally shaped wave relaxes back to a plane wave and no fold 
develops. Snapshots of a movie [31] with (a) t = 43, (b) t = 67, (c) t = 159, 
and (d) t = 400. The values of the feedback parameter are amin = 0-05 and 
ttmax = 0.15 and the control is switched on at t = 40. 


characterised by a continuing segmentation of waves and spreading spiral turbulence 
as shown in figure 2. For the same parameter values, patterns stop to segment after the 
feedback is switched on, giving rise to a persistent plane wave and two counter rotating 
spiral waves, see figure 12. The wave front of rotating patterns has a nonbinding 
positive curvature. According to the linear eikonal equation (13), it advances slower 
than a plane wave if the effective curvature coefficient z> is positive. Therefore, the 
plane wave has a tendency to annihilate rotating waves, finally leading to a solitary 
plane wave. 
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Figure 12. Curvature-dependent feedback suppresses spiral turbulence. After 
the feedback control is switched on at t = 256 (b), waves stop to break up, 
leaving behind a plane wave and a pair of counter rotating spiral waves (c) 
and finally a solitary plane wave (d). Same parameters deep in the regime of 
transversal instabilities as the corresponding uncontrolled time evolution in figure 
2. Snapshots of a movie [31] with (a) t = 5, (b) t = 301, (c) t = 672, and (d) 
t = 891. The values of the feedback parameter are amin = 0.04 and Umax = 0.08. 


4. Conclusions 

In this article, we present a feedback loop to induce, control, and suppress transversal 
instabilities of reaction-diffusion waves. The control signal is calculated from the 
local curvature of the isoconcentration line of the wave. We show that the curvature 
dependent control can amplify or quench small curvature perturbations in the wave 
shape. Simultaneously, the feedback allows to study a large variety of artificially 
produced wave patterns associated with transversal instabilities. Often these patterns 
are non-stationary and sensitively depending on small changes in the initial conditions 
characteristic for chaotic dynamics. 

Mathematically, the onset of transversal instabilities can be understood with the 
help of the linear eikonal equation which relates the wave velocity normal to an 
iso concentration line to its local curvature. The coefficient n in front of the curvature 
determines the stability of a flat wave. For positive values of z/, convex wave segments 
slow down while concave wave segments propagate at a higher velocity. Under these 
conditions a perturbed flat traveling wave recovers its flat shape. In the case of 
negative z/, a small positive curvature causes an increase of the wave velocity which 
in turn results in an increase of the local curvature. Now, a flat wave is unstable 
with respect to small curvature perturbations. The proposed feedback loop is able to 
change the sign of the coefficient z/. 

With experiments on chemical waves in the PBZR in mind, for realistic parameter 
values we show in numerical simulations with the Oregonator model that transversal 
instabilities of planar waves can be induced by the feedback. Right beyond the 
transversal instability of planar waves, we find nearly flat folded waves which are 
stationary in a comoving frame of reference. For weak feedback we observe small 
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ripple-shaped undulations traveling along the wave front. Upon increasing the 
feedback strength further, V-shaped wave patterns with spatio-temporal transversal 
modulations appear. These V-shaped waves travel at a velocity that depends on the 
opening angle but is considerably faster than that of the planar wave. Far away from 
the instability threshold, breakup of waves causes persistent annihilation and merging 
of excited domains, self-sustained rotatory motion and nucleation of rotating wave 
segments. Qualitatively, the emerging wave patterns correspond to those observed in 
numerical simulations with separated activator and inhibitor diffusivity [11]. 
Regarding chemical wave propagation in the PBZR, we emphasise that the feedback 
parameters of the control law are experimentally accessible. For appropriate BZ 
recipes the dependence of the wave velocity on the intensity of applied light should be 
strong enough to induce transversal wave instabilities. The isoconcentration line of 
the wave can be determined by 2d spectrophotometry with sufficient spatial resolution 
using the contrast between the oxidised and reduced form of the catalyst. We believe 
that the computation of the curvature by the Level Set Method as described in Sec. 2.4 
will work reliably for noisy experimental data, too. Because all chemical components 
share similarly shaped isoconcentration lines, the measurement of the concentration 
field of an arbitrary single chemical species is sufficient for setting up the control 
loop. Fine-tuning the feedback parameters allows to study the onset of transversal 
instabilities in dependence of the boundary conditions as e.g. the channel width L, as 
pointed out in chapter Sec. 2.2. 

In the opposite case, sufficiently strong feedback changes the sign of the effective 
curvature coefficient from negative to positive. Consequently, naturally occurring 
transversal wave instabilities leading to the breakup of waves are suppressed - the 
feedback stabilises planar waves and spiral waves. Spreading of spiral turbulence is 
inhibited due to the suppression of segmentation of waves. 

Reaction-diffusion waves describe, at least approximately, a huge variety of wave 
processes in biology. Our results are potentially applicable to deliberately induce or 
inhibit transversal wave instabilities and to control the emerging patterns under very 
general conditions. The essential condition for applicability is that the propagation 
velocity of the wave can be externally controlled over a sufficiently large range such 
that the curvature coefficient of the eikonal equation switches its sign. 

Moreover, we expect that curvature dependent feedback might have interesting 
applications in inter facial pattern formation in general. For example, this feedback 
mechanism could be the starting point for a control strategy aimed at the purposeful 
selection of patterns affected by interfacial instabilities as, e.g., alloys growing into an 
undercooled melt. 
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parameter 

value 

description 

/ 

1.4 

stoichiometric parameter 

Q 

0.002 

system parameter 

l/e 

49 

time scale separation 

1 /e 

4410 

time scale separation 

<l>o 

0.02 

background illumination 

Du 

1.0 

activator diffusion coefficient 

Du, 

1.2 

inhibitor diffusion coefficient 

V 

1.05 

curvature coefficient 

A 

0.68 

fourth order coefficient in the KS equation 

Cl 

-90.19 

slope of linear fit for velocity over 4> 

Co 

9.01 

constant of linear fit for velocity over 4> 

^norm 

1.2 

curvature normalisation 

At 

0.0001 

time step 

Ax,Ay 

0.05 

step width of spatial resolution 


Table Al. Parameter values used for numerical simulations of the modified 
Oregonator model. 


parameter 

value 

description 

a 

0.1 

excitation threshold 

kf 

2 

system parameter 

kg 

2 

system parameter 

a 

2.1 

ratio of diffusion coefficients 

s 

0.01 

system parameter 

e 

0.1575 

time scale separation 

Cl 

-8.75 

slope of linear fit for velocity over a 

Co 

2.23 

constant of linear fit for velocity over a 


Table A2. Parameter values used for numerical simulations of the piecewise 
linear FitzHugh-Nagumo model. 


Appendix B. Bonnet’s formula 

We prove the formula of Bonnet, i.e., we demonstrate that evaluating the curvature 
field defined by equation (27) at an isoconcentration line 7 yields the curvature of 7 . 
Let u = u{x^ y) be the map ^ M and 7 {y) = (7^, {y) , y)^ be the isoconcentration 
line 7 parametrised by It follows that u (7^, fy) ,y) = Uc = const, for 

all values of y. Therefore we can write (to shorten the notation, we write 

dx'^ix, y) (y), y)) 

-^u (7^ (y), y) = dxu (7^ (y), y) 7^ (y) + dyU (73, (y), y) = 0 , (B.l) 

cfi d 

{^x (y) ’ y) = ^ (y) ’ y) y'x (y) + 9yu (7® (y), y)) 


and 


































Control of transversal instabilities in reaetion-diffusion systems 


20 


= d:^u ( 7 ^ {y ), y) if (y) + d^u ( 7 ^ (y ), y) ( 7 ^ {y)f 
+ 2dyd^u ( 7 ^ {y ), y) 7 ^ (y) + ( 7 a; {v), v) 

= 0, (B.2) 

and generally -^u {-jx (y) ; y) = 0 with n G N, n > 0. The curvature field k, Eq. (27), 
expressed in Cartesian coordinates is 


K {x, y) = 


dlu {x, y) {dxU (x, y)f - 2dxU {x, y) dyU {x, y) dx,yU (x, y) 


({dxU {x, y)f + {dyU {x, y)f^ 


3/2 


dlu{x,y) (dyu(x,y)f 


(^{dxU {x, y)f + {dyU (x, y)f'^ 
Evaluating k at the isoconcentration line yields 


3/2 • 


(B.3) 


K iix (y), y) = 


dyU ( 7 a; (y), y) {dxU {jx (y), y)f + d^u {jx (y), y) {dyU {jx (y), y))" 


({dxU (7a; (y), y)f + {dyU {jx (y), 2/))^) 
25a,M (7a; (y), y) dyU (jx (y), y) dx,yU {-jx (y), y) 


X 3/2 


3/2 


(B.4) 


({dxU (7a; (y), y)f + (dyU {jx (y), y))^) 

Using Eq. (B.l), the denominator of Eq. (B.4) can be simplified as 

{dxU {jx (y), y)f + {dyU {jx (y), y)f = {dxU {jx (y), y)f + {-dxU (jx (y), y) ix {y)f 

= {dxu{^x{y),y)?{i + {ix{y)f)- (B.5) 

Similarly, using Eq. (B.l) and Eq. (B.2), the first term of the numerator of Eq. (B.4) 
can be rewritten in the form 

^7 ('1'® ( 2 /); y) (^^“ (7a; {y) , y)f = - {dxU {jx (y) , y)f (dxU (jx (y) , y) 7 " {x) (B.6) 

+dlu ( 7 a; (y), y) ( 7 ^, {x)f + 2dx,yU {^x (y), y) lx 

while the second term of the numerator of Eq. (B.4) can be cast as 

dlu {lx (y), y) {dyU (7a; (y), y)f = dlu {jx (y), y) {dxU {ix (y) , y)f (ix {x)f 
The last term of the numerator of Eq. (B.4) becomes 


(B.7) 


- ‘^dxU ( 7 a; (y), y) dyU {jx (y), y) dx,yU {jx (y), y) = 

2 {dxU ( 7 a; (y), y)f dx,yU {-/x (y), y) lx (y) ■ (B.8) 

All terms except the term proportional to if {x) in the numerator cancel. We are left 
with 

«( 7 . iv) ,y) = - - - , 3 / 2 > (B-9) 

(i + (7;(y))") 

which is exactly the curvature of a graph, see Eq. (15). 
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